MODULE VISCVAR
!
! Verwendungszweck:
!
!      Bereitstellung von Hilfsvariablen fuer die Viskositaet
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!       alle viscosen Hilfsgroessen sind zeitlich zentriert - wenn auch ohne 'Z' im Variablennamen geschrieben
!
!   divU        : Divergenz der Geschwindigkeit
!   divU_flag   : Schalter = 1 wenn divU < 0 , sonst 0
!   q_lin_r     : Parameter for linear viscosity                  \
!   q_quad_r    : Parameter for quadratic viscosity               |_  r-component
!   muQ_r       : scaling factor for radial viscosity             |   of viscosity
!   Uterm_r     : a term with velocities for the radial viscosity /
!   q_lin_p     : Parameter for linear viscosity                  \
!   q_quad_p    : Parameter for quadratic viscosity               |_  phi-component
!   muQ_p       : scaling factor for radial viscosity             |   of viscosity
!   Uterm_p     : a term with velocities for the radial viscosity /
!   q_lin_p_ene : \
!   alpha_p_ene : |-- separate variables to enable linear sheer viscosity in relaxation calculations
!   muQ_p_ene   : /
!
!-----------------------------------------------------------------------

      use config, only : rkind, np

      implicit none

      private
      public :: VISCO, q_lin_r, q_quad_r, alpha_r, q_lin_p, q_quad_p, alpha_p, &
      divU, divU_flag, muQ_r, Uterm_r, muQ_p, Uterm_p, &
      q_lin_p_ene, alpha_p_ene, muQ_p_ene


      real(kind=rkind) :: q_lin_r, q_quad_r, alpha_r(np), muQ_r(np), Uterm_r(np)
      real(kind=rkind) :: q_lin_p, q_quad_p, alpha_p(np), muQ_p(np), Uterm_p(np)
      real(kind=rkind) :: q_lin_p_ene, alpha_p_ene(np), muQ_p_ene(np)
      real(kind=rkind) :: divU(np), divU_flag(np)

CONTAINS





SUBROUTINE VISCO
!
! Verwendungszweck:
!
!     Berechnung von Hilfsgroessen fuer die Viskositaetsterme
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!   idx         : Schleifenindex
!

      use primvar, only : XZ, MR, MD, MUr, MUp
      use config,  only : np, BCflag
      use physco,  only : pi, z0, z1, z2, z4, z12, z23, grav
      use geomvar, only : S_volZ
      use matvar,  only : Pgas0Z
      use zvar,    only : H_pZ, rho0Z
      use RBvar,   only : M_star, R_in, R_out


      implicit none

      integer     :: idx

!-----------------------------------------------------------------------
!    alpha-model viscosity
!-----------------------------------------------------------------------

      alpha_r(:)     = q_lin_r     ! xxx man hat hier noch die Freiheit kompliziertere Modelle einzubauen, z.B. die Ionisation des Gases beschreibend
      alpha_p(:)     = q_lin_p
      alpha_p_ene(:) = q_lin_p_ene


!-----------------------------------------------------------------------
!    Hilfsvariablen der Viskositaet
!-----------------------------------------------------------------------
!
!  -> am Innenrand (idx=3) wird u_r geloest ... brauche viskose Hilfsvariable fuer idx = 2 & idx = 3
!  -> am Aussenrand (idx=np-2) wird u_r geloest ... brauche viskose Hilfsvariable fuer idx = np-2 & idx = np-3
!
!  -> in der innersten skalaren Zelle (idx=3) wird u_phi geloest ... brauche viskose Hilfsvariable fuer idx = 3 & idx = 4
!  -> in der aeussersten skalaren Zelle (idx=np-3) wird u_phi geloest ... brauche viskose Hilfsvariable fuer idx = np-2 & idx = np-3
!
!
!   Sinnvoll ausrechnen lassen sich die (skalaren) Viskositaetsvaribale nur von idx = 3 bis idx = np-3.
!   Fuer idx =2 und idx = np-2 gibt es zwar keine nummerischen Fehler, aber auch keine sinnvollen Werte,
!   letztendlich weil der viskose Druck in den Pseudozellen prinzipiell nicht definiert ist. In den
!   Pseudozellen wird daher der viskose Druck auch als P_visc = 0 annommen.
!
!   Uterm_p ist eine vektorielle Groesse, brauche ich fuer alle Gitterpunkte inkl. Randpunkte; d.h. von idx = 3 bis idx = np-2.
!

      !$OMP PARALLEL DO PRIVATE (idx)
      do idx = 2,np-2
         divU(idx) = z2*pi * ( XZ(MR,idx+1) * XZ(MUr,idx+1) - XZ(MR,idx) * XZ(MUr,idx) ) / ( S_volZ(idx) + BCflag(idx) )
      end do
      !$OMP END PARALLEL DO

      ! Flags fuer das Vorzeichen der Divergenz von u
      where ( divU(:) < z0 .and. BCflag(:) /= z1 ) ! not in pseudo-cells
         divU_flag(:) = z1
      elsewhere
         divU_flag(:) = z0
      end where


      !$OMP PARALLEL DO PRIVATE (idx)
      do idx = 2,np-2
         ! lineare Viskositaet mit Dicke der Scheibe ("alpha model")
         ! quadratische weiter mit R skaliert weil das fuer die Gitterweite relevant bleibt
         ! eigene Variable fuer Energiegleichung damit man Scherviskositaet zum relaxieren einschalten kann

         muQ_r(idx)     =  alpha_r(idx)     * H_pZ(idx) * XZ(MD,idx) * sqrt( Pgas0Z(idx) ) / sqrt( rho0Z(idx) ) &
                        - (q_quad_r * XZ(MR,idx))**2 * XZ(MD,idx) * divU(idx) * divU_flag(idx)

         muQ_p(idx)     =  alpha_p(idx)     * H_pZ(idx) * XZ(MD,idx) * sqrt( Pgas0Z(idx) ) / sqrt( rho0Z(idx) ) &
                        - (q_quad_p * XZ(MR,idx))**2 * XZ(MD,idx) * divU(idx) * divU_flag(idx)

         muQ_p_ene(idx) =  alpha_p_ene(idx) * H_pZ(idx) * XZ(MD,idx) * sqrt( Pgas0Z(idx) ) / sqrt( rho0Z(idx) ) &
                        - (q_quad_p * XZ(MR,idx))**2 * XZ(MD,idx) * divU(idx) * divU_flag(idx)


         ! in den Pseudo-Zellen idx = 2 und idx = np-2 haben die Radien XZ(MR,idx+1) und XZ(MR,idx) den gleichen Wert -> Schalter im Bruch + Uterm mit BCflag 0-setzen
         Uterm_r(idx) = (z1 - BCflag(idx)) * (   ( XZ(MUr,idx+1) - XZ(MUr,idx) ) / ( (z1 + BCflag(idx)) * XZ(MR,idx+1) - XZ(MR,idx) ) - &
                                                 ( XZ(MUr,idx)*XZ(MR,idx) + XZ(MUr,idx+1)*XZ(MR,idx+1) ) / ( XZ(MR,idx)**2 + XZ(MR,idx+1)**2 )   )

      end do
      !$OMP END PARALLEL DO

      !$OMP PARALLEL DO PRIVATE (idx)
      do idx = 3,np-2
         ! Kein viskoser Druck fuer die Randpunkte idx = 3 und idx = np-2 -> Uterm mit BCflag 0-setzen

         Uterm_p(idx) = (z1 - BCflag(idx) - BCflag(idx-1) ) * &
                        ( z2 * ( XZ(MUp,idx) - XZ(MUp,idx-1) ) / ( XZ(MR,idx+1) - XZ(MR,idx-1) ) - z12 * ( XZ(MUp,idx) + XZ(MUp,idx-1) ) / XZ(MR,idx) )

      end do
      !$OMP END PARALLEL DO

!-----------------------------------------------------------------------
!    Scherung in einer Kepler-Scheibe
!-----------------------------------------------------------------------
! Uterm = du/dr - u/r;  hier u = sqrt(GM/r); damit Uterm = -3/2 sqrt(GM/r^3)
!
! Fuer idx = 3 und idx = np-2 wird Uterm (und seine Ableitungen) auf 0 gesetzt, d.h. ich kann Uterm
! auf einen beliebigen _konstanten_ Wert setzen und die Ableitungen stimmen immer noch.

      Uterm_p(3)    = -z23 * sqrt(grav * M_star / R_in**3)
!     Uterm_p(np-2) = -z23 * sqrt(grav * M_star / R_out**3)
	  Uterm_p(np-2) = Z0

! kann man so eigentlich nicht machen, weil die Ableitungen nicht berücksichtigt sind, ich hoffe aber auf kleine Abweichungen um trotzdem eine Lösung zu finden
!      Uterm_p(np-2) = ( muQ_r(np-3) * ( XZ(MR,np-2)**2 + XZ(MR,np-3)**2 ) ) / ( muQ_r(np-2) * ( XZ(MR,np-1)**2 + XZ(MR,np-2)**2 ) )  * Uterm_r(np-3)
!      WRITE(*,*) 'Uterm_p am Außenrand = ', Uterm_p(np-2) 
!      WRITE(*,*) 'Uterm_p 2 am Außenrand = ', -z23 * sqrt(grav * M_star / R_out**3)



END SUBROUTINE VISCO

END MODULE VISCVAR
